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We consider non-interacting fermions on a lattice and give a general result for the 
reduced density matrices corresponding to parts of the system. This allows to calculate 
their spectra, which are essential in the DMRG method, by diagonalizing small matrices. 
We discuss these spectra and their typical features for various fermionic quantum chains 
and for the two-dimensional tight-binding model. 



I. INTRODUCTION 

Density matrices have found an interesting new application in recent years. In the density-matrix 
renormahzation group (DMRG) method @ || a quantum system is built from smaller parts and the idea 
is to work with basis functions in these parts which are optimal for the combined system. These are 
the eigenfunctions of the reduced density matrices which have the largest eigenvalues Wn- Obviously, 
the procedure will work well if the eigenvalue spectrum drops rapidly, such that a small number of 
functions practically exhausts the sum rule = 1. For quantum chains this is indeed the case. The 

numerical calculations show a roughly exponential decrease of the eigenvalues H,^. Of course, this raises 
the question whether such spectra can be obtained explicitely for some solvable models. For non-critical 
systems this is possible by using the relation Q between the density matrices of quantum chains and the 
corner transfer matrices (CTM's) of the corresponding two-dimensional classical problems. In this 
way, the spectra for the transverse Ising chain , the XXZ Heisenberg chain 0] and a chain of coupled 
oscillators |8j could be determined in the thermodynamic limit and compared with DMRG calculations. In 
all these cases, one finds simple analytic expressions and, apart from degeneracies, a strictly exponential 
behaviour. This does not hold for the chiral three-state Potts chain |^ or for non-integrable models 
but qualitatively the spectra are similar. 
Given the importance of fermionic systems in general and also for DMRG applications, one would of 
course like to have results for this case, too. The transverse Ising chain can be viewed as a fermionic model, 
but the CTM approach does not make use of this and is limited to large non-critical systems. Therefore 
an alternative approach is necessary by which one can treat solvable fermion systems of arbitrary size. 
In the present communication we show how this can be done. The systems which we consider are non- 
interacting, such that the Hamiltonian can be diagonalized by a Bogoliubov transformation. Using an 
explicit form of the state in question (usually the ground state), we show that arbitrary reduced density 
matrices can be calculated exactly and have the general form exp^—Ti.). The operator Ti describes a 
collection of new, non- interacting fermions with single-particle eigenvalues Si. Apart from the different 
statistics, this is the same situation as for coupled oscillators p|,p^. The which determine the properties 
of the spectrum, follow from the eigenvalues of an M x M matrix, where M is the number of sites in 
the chosen subsystem. In general, they have to be calculated numerically. One should stress that the 
dimensionality of the system plays no essential role. 

In the following Section ^ we sketch the general method of computation which uses coherent fermionic 



states for calculating the necessary partial traces. In Section [II, we apply it to the transverse Ising chain 



and discuss the resulting spectra for a number of situations, including the critical case, the first excited 



state and related row transfer matrices. Section IV deals briefly with another one-dimensional problem. 



namely the spin one-half XY-chain in a field. This is interesting because it has a disorder point where 
the spectrum collapses. In Section ^ we turn to the physically most important case of a tight-binding 
model which we dicuss in two dimensions. We present spectra for systems of various sizes and shapes, 
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as well as truncation errors showing the difficulties in this case. Section finally, contains a summary 
and some additional remarks. Some technical details can be found in the Appendix. 



II. METHOD 



We consider Hamiltonians which are quadratic in Fermi operators and thus have the general form 

i 1 

H = J2 + ^{4b,,c] + h.c.)] (1) 

where the Ci's and c|'s are Fermi annihilation and creation operators. Because of the Hermiticity of H, 
the matrix A is Hermitian and B is antisymmetric. In the following we consider only real matrices. One 
can diagonalize H through the canonical transformation [l^ 

Vk = ^(fffeiQ + hkicl) (2) 

i 

which leads to 

H = Akijlijk + constant (3) 

k 

The quantities are the eigenvalues of the matrices {A — B){A + B) and {A + B){A — B), the 
corresponding eigenvectors being (pki = gki + hki and ipki = 9ki — /ifei, respectively. 

Consider now the ground state | $o > of the Hamiltonian (nl) for an even number of sites L. Due to the 
structure of H, it is a superposition of configurations with either an even or an odd number of fermions. 
This suggests to write it (for the even case) in the form 



|$o>= C exp{i^G.,clc]}|0> (4) 
where | > is the vacuum of the Ci, i.e. 
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I >= (5) 

Such an exponential form is known from superconductivity, where the BCS wave function (in momentum 
space) can be written in this way p^ . 

One obtains Gij by applying the Fermi operators rjk to the ground state 

r^fe I $0 > = for aU k (6) 

which leads to (see Appendix) 

^ 9kmGmn + h^n = for all fc, n (7) 

Thus G relates the two matrices g and h of the transformation (|^). Using (|^), one obtains the total 
density matrix po =| >< 'I'o | explicitly in an exponential form 

po = exp (i G^^clc]) I X I exp (-^ ^ G.,c,c,) (8) 

ij ij 
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One now divides the total system in two parts (system and environment in the DMRG terminology) 
and looks for the reduced density matrix in part 1. This is obtained by taking the trace over part 2. 

Pi = Tr2 (po) (9) 
In order to calculate pi , one uses the fermionic coherent states defined by 

c. >-e. > (lo) 

Such states can be built from the vacuum with operators ci and Grassmann variables 

|^i...ei>=exp(-^^,cl)|0> (11) 

i 

Using this, one can write the trace of an operator O as 

TtO= [U dCd^^e- S„ < _^ I o I ^ > (12) 

a 

After forming a general matrix element of po with such states and taking the trace over the environment 
with ([l2|), one obtains, if part 1 consists of M sites 



< 6 • • • a/ I Pi I • • • Cm > 

= \C\' n < 6 • ■ • - Cm+1 • ■ • - a I po I • ■ • ■■<l> (13) 

Inserting (^) leads to an integrand which contains only quadratic forms of Grassmann variables in the 
exponents. The integration can then be carried out by rotating and displacing the variables as for a 
Gaussian integral with complex numbers. This gives 

< Ci ■ ■ • 6/ I Pi I • ■ • Cm > 

= |Cp exp {J2 exp A.ere^) exp (^ -a.,C.'C;) ; ^,j<M (14) 

ij ij ij 

The M X AI matrices a and (3 appearing here are defined as follows. One divides G into four submatrices 
a^^,a^^,a'^^ and a^^, according to whether the sites i,j belong to part 1 or part 2. In terms of these 

2 a = a" + ca^^c^ 
(3 = cJ (15) 

where c = a^^(l — a^^)^^ and denotes its transpose. As shown in the Appendix one can reconstruct 
the operator form of pi from the matrix elements ( |l4|) . This gives 

pi = |Cp exp a^jcjc]) exp (^(ln/3)ijc|cj) exp (^ -a^jdCj) ; i,j <M (16) 

ij ij ij 

Finally, since the Fermi operators appear quadratic in the exponents, pi can be diagonalized with a 
Bogoliubov transformation as in (||). As a result, 

M 

pi =X exp (-^ £,///,) (17) 
1=1 
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with new Fermi operators //, /; and K = |Cp. The single-particle eigenvalues ei follow from the matrices 
a, (3 according to Eqn. ( |34| ) of the Appendix. The normalization factor K is fixed by the sum rule 
Tr(pi) = 1. In this way, one can calculate the density-matrix spectra numerically for an arbitary part of 
a finite system with Hamiltonian (Fl). 




FIG. 1. Single-particle eigenvalues ei for one half of a transverse Ising chain, arranged in ascending order. The 
system is in the groundstate, L = 20 and A < 1. 
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FIG. 2. Density-matrix eigenvalues ui„ , arranged in decreasing order, obtained from the e; in Fig. |l] and for 
the same parameters. 



III. TRANSVERSE ISING CHAIN 



by 



As a first example,we consider in this section the transverse Ising chain with open boundaries described 



L-l 



(18) 
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where the cr" are PauU spin matrices and the transverse field has been set equal to one. In the thermo- 
dynamic limit, this system has a quantum critical point at A = 1 and long-range order in cr^ for A > 1. 
In terms of spinless fermions H reads 

L L-l 

H = -2^(clc. - 1/2) - A ^(4 - c.)(4+i -I- c.+i) (19) 

1=1 i=l 

and thus has the form (|l|). In the following we discuss the reduced density matrix pi for one half of the 
chain, i.e. M = L/2. 

We first consider the ground state. In figure the single-particle eigenvalues £i are plotted for L = 20 
and different coupling constants A. For A = 0.1 they all lie on a straight line, which corresponds to the 
situation one finds in the thermodynamic limit. This is what one expects since the correlation length is 
much less than L and hence boundary effects should be small. One can also check that the values are 
exactly those obtained analytically via corner transfer matrices [Q. It seems to be difficult, however, to 
derive these results directly from our equations. For larger coupling, A = 0.5, only the first ei follow 
a linear law, then the curve bends upwards. This is similar to the behaviour one finds for finite-size 
corner transfer matrices |]l6| , although the geometry there is different. At the same time, the initial slope 
decreases. Finally, at the critical point, the whole graph is curved. In the ordered region (not shown), a 
linear regime develops again. 

From the si one obtains the actual eigenvalues Wn of pi by specifying the occupation numbers in 
(p^). The resulting spectra are shown in figure ^ in a semi- logarithmic plot. Note that not all w„ are 
shown, however they are correctly normalized to one. Similar results, but for a smaller number of w„, 
were obtained in via DMRG calculations. Due to the relatively large values of the ei there is a rather 
rapid decay (note the vertical scale) so that the system can be treated very well by DA/IRG [ pT|JT^ . This 
holds even at the critical point, where the decay is slowest. 

The situation there is presented in more detail in the next figures. Figure ^ shows the e-spectra 
for various sizes of the system. As L increases, the number of e increases, the curves become ffatter, 
but the curvature remains. There is no sign of a linear region related to conformal invariance on this 
scale (compare [|l6|). The w„ spectra are plotted in figure ^ Because of the form of the e, there are 
few degeneracies and the curves have the typical, relatively smooth shape found also for other critical 
systems [|,|. The finite -size effects show up essentially in the tails. 
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FIG. 4. Density-matrix eigenvalues w„ for transverse Ising chains at the critical point obtained from the ei in 
Fig. I 

So far, we have treated the ground state, but one can also determine the density matrices for the first 
excited state | <i>i >. This state contains an odd number of fcrmions. To apply the formalism here, one 
can perform a particle-hole transformation at one site, e.g. c\ ^ ci. Then | $1 > appears in the even 
subspace and can be written in the form (^. With the help of the relations 

ryt I $1 > - 

77fe I $1 > = for /c > 2 (20) 

one can then derive the corresponding equation for the matrix Gij. In this way, the single-particle 
eigenvalues ei shown in figure ^ were obtained. In contrast to the case of the ground state, the first 
eigenvalue is zero here. This refiects the fact that, in the original representation, the fermion number is 
odd, while the number of sites is even. The other eigenvalues are very similar to those for the ground 
state. In particular, one has a linear spectrum away from A = 1 and a curved one at the critical point. 
The vanishing ei causes all eigenvalues w„ of pi to be at least doubly degenerate. 




FIG. 5. Single-particle eigenvalues si for the first excited state of a transverse Ising chain for L = 12 and four 
values of A. 

Finally, the closely related problem of the row-to-row transfer matrices for the two-dimensional Ising 
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model can be studied in the same way. For a square lattice with couplings Ki {K2) in the vertical 
(horizontal) direction one can consider two symmetrized versions, namely 



V ^V2^''^ViV2^'^] W = Vi''^V2Vi''^ (21) 
where Vi (V2) contain the vertical (horizontal) bonds. Both represent ferniionic quantum chains and can 



be diagonalized also for open boundaries 20|. For the thermodynamics, one needs the eigenvector 
with maximal eigenvalue. DMRG calculations using the operator V have already been done Ip^- The 
spectrum of the £; in the isotropic case Ki = K2 is very similar to that found above in figurcll]. This 
also holds for the magnitude of the e; and the problem can therefore be treated equally well by DMRG. 
For W the e-spectrum is strictly linear at the lower end and described by a formula containing elliptic 
integrals as in , while for V the values are somewhat smaller and there is a deviation from linearity for 
the first ei. This reflects the difference in the representation of pi via CTM's in the two cases. 



IV. XY-SPIN CHAIN 



In this section we consider briefly the spin one-half quantum chain described by the Hamiltonian 



L-1 



H = - J/2 ^[(1 + 7)a,ra,r+i + (1 - j)afaf^, + h{a^ + a^,)] (22) 

i=l 

which reads in terms of fermions 

L-1 

H = -jJ2{ich+i + 744+1 + h.c.) + h{cU + 4+iC^+i - 1)] (23) 

i=l 

Although similar to the transverse Ising chain, this system has a special feature. For 

7^ + 1 (24) 

the ground state simplifies and also becomes two-fold degenerate. In the spin language, one has two 
simple product states |^ . Moreover, the behaviour of correlation functions changes from monotonic to 
oscillatory and thus (24) represents a "disorder line" |Q. On this line, H describes also a stochastic 
reaction-diffusion model |25| equivalent to Glauber's kinetic spin model. 

The appearance of a simple ground state can be observed in the density-matrix spectrum and has 
already been seen in DMRG calculations for certain other models (see section 3.1 in [|j). For the XY- 
chain, it can be investigated very well in the fermionic approach. 

In figure |^ we show the lowest e;-values as a function of the parameter h for fixed 7 = 1/2. The 
disorder point according to ( p4| ) is then at ho = 0.866. One can see that, coming from larger values of ft., 
all £1 except the lowest one diverge as one approaches h^. For h < ho they become finite again. In this 
region, however, one has to work in another subspace since at ho the lowest fermionic eigenvalue Aq in 
H crosses zero, which leads to the degeneracy of the ground state. This can be done as for the excited 



state in Section III . Then one finds the curves in the figure. As a check we also performed direct DMRG 
calculations and found complete agreement (dots). Such crossings appear repeatedly as one reduces h 
further. The next one (for the chosen L) takes place at ft, = 0.78. However, as seen from the figure, the 
higher ei show no effects at this point, indicating that the ground state of H does not simplify there. At 
/iQ, the divergence of the £; for I > 2 together with the value ei = lead to the density-matrix eigenvalues 
wi — — 1/2, while all other w„ are zero, i.e. the spectrum collapses at this point. This effect could 
be a tool in the search for simple ground states by DMRG. 
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FIG. 6. The four lowest single-particle eigenvalues e for an XY spin chain in a field h. The anisotropy is 7 = 0.5, 
the length L = 8. Lines result from the analytical method, solid circles from a DMRG calculation. 



V. TWO-DIMENSIONAL TIGHT-BINDING MODEL 



As the last, but most important example we consider a tight-binding model with open boundaries 
described by 

^ = - E(^i^j+4^0 (25) 

<ij> 

where the brackets < i,j > denote nearest-neighbour sites. This model is critical and solvable in all 
dimensions. We treat it here for the case of a square lattice and we assume that the system also has the 
shape of a square with L = N'^ sites where N is even. This problem has has served as a DMRG test case 
some time ago p6[ . 

The ground state here is different from that in the previous sections. Because H only contains hopping, 
_B = in Equ.(|l|), the fermion number is fixed and | $o > does not have the form (^. However, one 
can perform a particle- hole transformation on L/2 sites, for example on every second one, by which 
the Hamiltonian acquires pair creation and annihilation terms (B ^ 0). Then | $0 >: which originally 
contains L/2 particles, becomes a superposition of terms with particle numbers ranging from to L and 
can again be written in the form . In the same way, an arbitrary n-particle eigenstate of H could be 
handled by exchanging particles and holes at n sites. The density-matrix spectrum is not affected by 
such local transformations. 

To carry out the calculation, one makes the problem formally one-dimensional by numbering the sites 
from 1 to L in such a way that the desired partition into two parts arises naturally. For example, a 
meander-like numbering as in p6[ | permits a division of the square into two halves. 

In figure 0, the single-particle eigenvalues ei for such a half-system and three different sizes are shown. 
One notices two features which are in contrast to the one-dimensional results: a "foot" of low- lying ei 
and a much smaller slope of the curves (note the scales). Both are strongly size-dependent. The number 
of ei in the foot is equal to N, which indicates that these states are closely connected with the interface 
between system and environment. Figure ^ shows the first 2000 eigenvalues w„ which result. Due to the 
small El, they decrease very slowly and the situation worsens as the system is enlarged. The tails of the 
curves can be described qualitatively by ln(iy„) ~ — In^(n) as in | 
up even more in the truncation error /„, which is defined as the sum of all w's beyond n. This quantity 
is given in the inset of the figure. With n = 2000 it is approximately 5 x 10^^, 5 x 10^^ and 10~^, 



11 13]. The effect of these tails shows 
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respectively. Thus the situation is not only much worse than for one-dimensional systems, but also worse 
than for the two-dimensional system with a gap discussed in |l2j. Standard DMRG calculations using, 
say, 2000 states would be limited to sizes below 12 x 12, and even then the accuracy would be much less 
than one is used to in quantum chains. 
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FIG. 7. Single-particle eigenvalues e; for two-dimensional tight-binding models of different sizes. The ei are for 
one half of the system. 
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FIG. 8. Density-matrix eigenvalues ui„ of two-dimensional tight-binding models, obtained from the ei in Fig. 
The inset shows the truncation error (see text). 



One can also calculate the density-matrix spectra for other shapes of the selected subsystem. As an 
illustration, we show in figure ^results for one quarter of a quadratic system (for example the upper right 
one). Note that the sizes indicated there refer to the whole system. One sees again some small eigenvalues, 
but fewer than for the half-system, while there are further higher-lying plateaus and additional short steps. 
Obviously this reflects the particular interface with a corner. For the 10 x 10 system, for example, the 
two lowest plateaus contain 9 states which is just the number of sites along the interface. The eigenvalues 
Wn are plotted in the inset of the figure. They are similar to those for the half-system but some more 
steps persist for small n. In the same way, one can investigate cases where one cuts the square diagonally 
at various positions. Such partitions appear in a recent new DMRG algorithm |2^ . The general features 
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of the spectra, however, do not change. 

Finally, let us mention that one can also include spin in H and thereby treat the Hubbard model in 
the U = hmit. Then the operators // in pi acquire a spin label, too, and all single-particle levels 
become doubly degenerate. This makes the tails of the ii;„-curves even flatter than in the spinless case. 
However, the curves are also pulled down by smaller normalization factors which leads to a faster initial 
decay. For a 20 x 20 lattice, the spectrum of the first 3000 states is, on the whole, rather close to that 
shown in figure |[ 
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FIG. 9. Single-particle eigenvalues ei of two-dimensional tight-binding models. The £i are for a quarter of the 
system, the Wn obtained from them are plotted in the inset. 



VI. CONCLUSION 



We have studied the reduced density matrices for non-interacting fermions on a lattice. The key 
ingredient for the calculation was a simple representation of the (ground) state. This led rather directly to 
the exponential Boltzmann-like form of the density matrices. The only really numerical step involved was 
the calculation of the single-fermion eigenvalues appearing in the exponent. With these, we discussed a 
number of cases in one and two dimensions with characteristic differences. We focused on the eigenvalues, 
but one can also investigate the single- fermion eigenfunctions. One then sees that they are concentrated 
near the interface between the two parts of the system. This explains the decisive role of the connectivity 
for the spectra. 

One should mention that fermionic density matrices have been studied before, e.g. in quantum chem- 
istry mjl^. However, in this case the systems are continuous and the Hilbert space is infinite. Then 
already the single-particle density matrices have inifinitely many eigenstates Our systems are dis- 
crete, but we are interested in density matrices for arbitrarily large subsystems. These are non-trivial 
even for non-interacting fermions. From the experience with other models, one can expect that the results 
are roughly representative also for more complicated systems. 

For this reason, the two-dimensional case is particularly important. With our formulae, we could treat 
the tight-binding model for arbitrary partitions of the system. This allows to make much more detailed 
statements than a previous, purely numerical investigation of this system p6| . In particular, one can see 
the very slow decay of the spectra and of the truncation errors directly. Basically, it is connected with 
the existence of long boundaries between the two parts of the system. In the current DMRG procedures, 
these appear necessarily at some point of the calculation. Therefore it is not yet clear whether a recent 
new algorithm [p7| can really overcome this problem. 
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APPENDIX 



Here we hst some details concerning the steps in section (|l|) 
(A) To derive (|^), one writes the Eqn. explicitly as 



J2igknCn + hknci)e^ | (26) 



where F = 1/2^- G^jcjct Using the relation 



d 



[c„en = fxe^ (27) 
del 

the exponential factor can be brought to the left 

e^^{^5fcmG™„ + /ifen}c], I 0>=0 (28) 

n m 

Since this must hold for all k, the only possibility is that the term in the bracket vanishes which gives 
the desired result. 

(B) The explicit form of the integrand in (|l^) is 

exp{-Cf 6 + - ^'^a''^) ^ {^fa^H2 + ^'^a'^CD + imfo^^'Ci " ^^""^O} (29) 

where ^1,^1(^2 '^2) are vectors composed of the variables of part 1 (part 2), respectively. Using the 
notation ^ = (^2,^2)1 this can be rewritten as 

exp{-^t5^ + ^t^ + ^t^ + ^} (30) 

where B is a 2{L — M) x 2(L — M) matrix containing a^^, C,,rj are both 2{L — M) dimensional vectors 
constructed from a^"^ , a?^ , 5,1 and and K is the last term in (^9|). ( ^o|) is an explicit Gaussian form 
which can be integrated whereby (|lj) is obtained. 

(C) To derive the operator form for pi from Eqn. (|l4|), one first diagonalizes the matrix (3. This 
transforms (p3) into a similar form with modified matrix a. Using the relations 



del I 



c,c, 1^^,^; > le-e; > (31) 

one can replace with cjc] and ^-^j with acj in the left and right exponentials. The cross terms 
e-^'5i 4i , where \i is one of the eigenvalues of /3, can be rewritten with the relation 

<6l/(c!,Q)|e^>=e«*«'/(C,eD (32) 
In our case the left-hand side equals e'^'^^'>* = 1 -}- A,;^*^- so that 
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/(cl,c.) = (l + (A.-l)4c.) 

= e''^^'"^'^* (33) 

Transforming back to the original representation leads to 

(D) The operator pi in ( |l6| ) can be diagonalized by calculating the Heisenberg operators piCjPi^ and 
as in [Q. Due to the form of pi, they are linear combinations of the c and c^. Inserting the 
Bogoliubov transformation and following ||2^ one finds that the eigenvalues e; can be obtained from the 
equation 

(/3 + (3-^ + p-^a - ap-^ - aP'^a) xi = 2 coshe, xi (34) 

Typically, the matix /3 has elements which vary exponentially over a large range. This limits the size of 
the systems for which one can use Eqn. ( |3^ ) in actual numerical calculations. 
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